Übungsblatt 2
Aufgabe 1
Der Datensatz no2basel.rda enthält Angaben über die Menge an Stickstoffdioxid in der Basler Luft. Als Prädiktoren kommen der Tag der Messung, die Temperatur und die Windstärke zum Einsatz.
Teilaufgabe a)
Passen sie ein geeignetes, multiples lineares Regressionsmodell an. Erwägen sie durch Analyse von Residuenplots, ob eine Transformation der Zielgrösse notwendig erscheint. Prüfen sie dann die partiellen Residuenplots. Welche Auffälligkeiten bestehen bei den Prädiktoren?
Call:
lm(formula = NO2 ~ ., data = no2basel)
Residuals:
Min 1Q Median 3Q Max
-27.805 -12.779 2.763 8.251 34.099
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 97.6304 13.4975 7.233 1.79e-07 ***
Tag -0.1349 0.6742 -0.200 0.84305
Temp 3.5446 1.1660 3.040 0.00564 **
Wind -9.9637 3.0395 -3.278 0.00318 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 16.05 on 24 degrees of freedom
Multiple R-squared: 0.6321, Adjusted R-squared: 0.5861
F-statistic: 13.75 on 3 and 24 DF, p-value: 2.012e-05
Call:
lm(formula = NO2 ~ Temp + log(Wind) + Tag, data = no2basel)
Residuals:
Min 1Q Median 3Q Max
-26.612 -12.223 1.964 7.313 33.777
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 87.8353 11.6486 7.540 8.85e-08 ***
Temp 3.5038 1.1584 3.025 0.00585 **
log(Wind) -18.0835 5.4245 -3.334 0.00277 **
Tag -0.1638 0.6723 -0.244 0.80956
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.97 on 24 degrees of freedom
Multiple R-squared: 0.636, Adjusted R-squared: 0.5905
F-statistic: 13.98 on 3 and 24 DF, p-value: 1.778e-05
Call:
lm(formula = log(NO2) ~ ., data = no2basel)
Residuals:
Min 1Q Median 3Q Max
-0.38837 -0.13873 0.05543 0.11652 0.46320
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.624286 0.192467 24.026 < 2e-16 ***
Tag -0.003261 0.009614 -0.339 0.73737
Temp 0.047775 0.016627 2.873 0.00837 **
Wind -0.151553 0.043342 -3.497 0.00186 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2289 on 24 degrees of freedom
Multiple R-squared: 0.611, Adjusted R-squared: 0.5624
F-statistic: 12.57 on 3 and 24 DF, p-value: 3.87e-05
Call:
lm(formula = log(NO2) ~ Temp + log(Wind) + Tag, data = no2basel)
Residuals:
Min 1Q Median 3Q Max
-0.53084 -0.10450 0.05867 0.12760 0.44721
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.456994 0.171170 26.038 < 2e-16 ***
Temp 0.046769 0.017022 2.748 0.01121 *
log(Wind) -0.257980 0.079711 -3.236 0.00352 **
Tag -0.003190 0.009879 -0.323 0.74956
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2346 on 24 degrees of freedom
Multiple R-squared: 0.5913, Adjusted R-squared: 0.5402
F-statistic: 11.57 on 3 and 24 DF, p-value: 6.914e-05
Teilaufgabe b)
Modellieren sie den NO2-Wert nun mit einem GAM.
Entscheiden sie, ob eine Log-Transformation der Zielgrösse nötig ist.
Prüfen sie die Form der flexiblen Komponenten. Scheint deren Einfluss realistisch, bzw. mit einem physikalischen Verständnis des Problems vereinbar? Greifen sie wo nötig ein und verändern sie die Glättungsparameter.
Führen sie ebenfalls (mit einer Methode nach Wahl) eine Variablenselektion aus, um zu entscheiden, ob gewisse Terme aus dem Modell eliminiert werden können, bzw. eliminiert werden sollen.
Entscheiden sie sich dann für eine bestes GAM.
Family: gaussian
Link function: identity
Formula:
NO2 ~ s(Temp) + s(Wind) + s(Tag)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 73.357 1.297 56.58 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Temp) 8.200 8.763 14.90 1.86e-05 ***
s(Wind) 4.782 5.732 13.71 7.31e-05 ***
s(Tag) 1.000 1.000 20.84 0.000534 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.924 Deviance explained = 96.4%
GCV = 101.25 Scale est. = 47.073 n = 28
Method: GCV Optimizer: magic
Smoothing parameter selection converged after 18 iterations.
The RMS GCV score gradient at convergence was 7.202671e-06 .
The Hessian was positive definite.
Model rank = 28 / 28
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(Temp) 9.00 8.20 1.44 0.96
s(Wind) 9.00 4.78 1.13 0.67
s(Tag) 9.00 1.00 1.12 0.66
Family: gaussian
Link function: identity
Formula:
log(NO2) ~ s(Temp) + s(Wind) + s(Tag)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.23867 0.02052 206.6 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Temp) 8.176 8.784 2.575 0.077769 .
s(Wind) 1.000 1.000 18.187 0.000925 ***
s(Tag) 5.107 6.099 5.369 0.005191 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.902 Deviance explained = 95.4%
GCV = 0.025961 Scale est. = 0.011791 n = 28
Method: GCV Optimizer: magic
Smoothing parameter selection converged after 23 iterations.
The RMS GCV score gradient at convergence was 1.277338e-07 .
The Hessian was positive definite.
Model rank = 28 / 28
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(Temp) 9.00 8.18 1.38 0.98
s(Wind) 9.00 1.00 1.09 0.60
s(Tag) 9.00 5.11 1.03 0.47| df | AIC | |
|---|---|---|
| fit.gam | 15.98226 | 197.82772 |
| fit.gam1 | 16.28294 | -34.40441 |
Achtung: die Modellvergleichsparameter GCV, AIC oder auch das Adjusted R-Squared helfen uns hier wegen der Transformation der Zielgrösse nicht weiter (d.h. die Werte sind nicht vergleichbar!).
i. Aber die Residuenplots sprechen dafür. Insgesamt ist keine systematischer Fehler erkennbar.
ii. Die gefitteten Kurven, vorallem bei der Temperatur zeigen einige Peaks, sind diese realistisch? oder ist dies overfitting? Im nächsten Schritt wird etwas mehr Glättung angepasst.
library(mgcv)
par(mfrow = c(2, 2))
fit.gam2 <- gam(NO2 ~ s(Temp) + s(Wind) + s(Tag), data = no2basel, gamma = 1.266)
plot(fit.gam2, residuals = TRUE, shade = TRUE, pch = 20, cex = 0.5)
par(mfrow = c(2, 2))
gam.check(fit.gam2, pch = 20, rep = 100)
Method: GCV Optimizer: magic
Smoothing parameter selection converged after 25 iterations.
The RMS GCV score gradient at convergence was 1.422634e-05 .
The Hessian was positive definite.
Model rank = 28 / 28
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(Temp) 9.00 1.00 0.69 0.03 *
s(Wind) 9.00 1.38 0.93 0.28
s(Tag) 9.00 4.67 0.77 0.07 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit.gam3 <- gam(log(NO2) ~ s(Temp) + s(Wind) + s(Tag), data = no2basel, gamma = 1.266)
plot(fit.gam3, residuals = TRUE, shade = TRUE, pch = 20, cex = 0.5)
par(mfrow = c(2, 2))
gam.check(fit.gam3, pch = 20, rep = 100)
Method: GCV Optimizer: magic
Smoothing parameter selection converged after 16 iterations.
The RMS GCV score gradient at convergence was 5.03779e-08 .
The Hessian was positive definite.
Model rank = 28 / 28
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(Temp) 9.00 1.00 0.66 0.020 *
s(Wind) 9.00 1.00 0.93 0.250
s(Tag) 9.00 4.57 0.80 0.085 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ii. Mit gamma von 1.266 werden die Fits für Temperatur und Wind linear. Die Variable Tag dagegen erhält einen polynom. Das Bild scheint physikalisch realistischer zu sein .
iii. Im nächsten Schritt wird eine Variablenselektion durchgeführt.
Method: GCV Optimizer: magic
Smoothing parameter selection converged after 109 iterations.
The RMS GCV score gradient at convergence was 5.805592e-05 .
The Hessian was not positive definite.
Model rank = 28 / 28
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(Temp) 9.00e+00 4.50e-07 0.65 0.03 *
s(Wind) 9.00e+00 9.47e-01 0.94 0.28
s(Tag) 9.00e+00 4.68e+00 0.77 0.09 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Method: GCV Optimizer: magic
Smoothing parameter selection converged after 97 iterations.
The RMS GCV score gradient at convergence was 4.744154e-06 .
The Hessian was positive definite.
Model rank = 28 / 28
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(Temp) 9.00e+00 3.84e-10 0.69 0.045 *
s(Wind) 9.00e+00 1.39e+00 0.96 0.315
s(Tag) 9.00e+00 4.82e+00 0.72 0.025 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Teilaufgabe c)
Nutzen sie die Erkenntnisse aus b), aber ziehen sie nun auch noch Interaktionsterme in Erwägung. Ziehen sie ein abschliessendes Fazit, ob Interaktionen nötig sind.
Teilaufgabe d)
Welches der bisherigen Modelle eignet sich am besten? Verwenden sie Modellvergleichskriterien und statistische Tests, um zu einer Entscheidung zu kommen.
Teilaufgabe e)
Erzeugen sie eine erwartungstreue Vorhersage mit einem 95%-Vertrauensintervall für Tag=7, Temperatur=3 und Wind=4.
Aufgabe 2
In dieser Aufgabe kommt der kleine (simulierte) Lohndatensatz in lohn.rda zur Anwendung, der auf den Slides bereits vorgestellt wurde. Enthalten ist der Monatslohn von 200 Personen, deren Alter sowie das Geschlecht (m/f). Wir studieren in dieser Aufgabe die unterschiedlichen Fits, die wir von einem flexiblen Modell mit/ohne Interaktionsterm erhalten.
load("lohn.rda")Teilaufgabe a)
Passen sie das additive Modell lohn ~ s(alter) + geschlecht an. Plotten sie dann die gefitteten Lohn-Kurven für Männer und Frauen in einem Scatterplot mit den Daten ein (Hinweis: mit predict() arbeiten!). Erzeugen sie auch einen Tukey-Anscombe-Plot um zu prüfen, ob das Modell passt. Welche Aussage ist aus dem Summary über die beiden Kurven möglich?
Teilaufgabe b)
Verwenden sie nun das Interaktions-Modell lohn ~ s(alter,by=geschlecht) + geschlecht. Plotten sie erneut die gefitteten Lohn-Kurven für Männer und Frauen in den Scatterplot ein. Was sagt ihnen das Summary hier?
Teilaufgabe c)
Prüfen sie nun, was lohn ~ s(alter, by=geschlecht) ergibt, sowohl via grafische Darstellung wie durch Studium des Summary-Outputs. Was sind die Erkenntnisse daraus?
Teilaufgabe d)
Führen sie mit der Funktion gam.check() eine Residuenanalyse für das Interaktions-Modell aus b) durch. Was sind ihre Erkenntnisse daraus und wie sollte das Modell modifiziert werden? Rechnen sie dann dieses modifizierte Modell sowohl mit als auch ohne Interaktionsterm und stellen sie die Fits in den Scatterplots auf der Originalskala (d.h. Lohn vs. Alter) dar.
Aufgabe 3
Wir verwenden in dieser Aufgabe einen Datensatz, wo es um die Modellierung von Herzkrankheiten bzw. Bluthochdruck gibt. Sie finden den Datensatz als heart.rda. Er hat 462 Beobachtungen und 10 Variablen.
| Variable | Beschreibung |
|---|---|
sbp |
Systolic blood pressure. |
tobacco |
Cumulative tobacco (kg). |
ldl |
Low-density lipoprotein cholesterol. |
adiposity |
Adiposity. |
famhist |
Family history of heart disease (Present, Absent). |
typea |
Type-A behavior. |
obesity |
Obesity. |
alcohol |
Current alcohol consumption. |
age |
Age at onset. |
chd |
Response, coronary heart disease. |
Die Variable sbp ist die Zielgrösse, die anderen Variablen stehen als Prädiktoren zur Verfügung. Folgende Aufgaben sind zu lösen:
Teilaufgabe a)
Verwenden sie zuerst nur die beiden Variablen obesity und age als Prädiktoren. Finden sie ein geeignetes Modell zur Beschreibung der Zielgrösse. Erwägen sie dabei sowohl lineare wie auch additive Modelle, allenfalls mit der Verwendung von Interaktionstermen. Vergleichen sie diese mit den Modellvergleichsparametern sowie auch mit statistischem Testen.
Teilaufgabe b)
Verwenden sie nun alle zur Verfügung stehenden Prädiktoren, reduzieren sie das Modell jedoch auf die wesentlichen Prädiktoren. Verzichten sie dieses Mal zwecks einfacherer Interpretierbarkeit auf Interaktionsterme, erwägen sie aber immer noch ein lineares Modell und ein GAM.